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1 Introduction 

Accurate numerical simulations of complex multiscale compressible viscous 
flows, especially high speed turbulence combustion and acoustics, demand 
high order schemes with adaptive numerical dissipation controls. Standard 
high resolution shock-capturing methods are too dissipative to capture the 
small scales and/or long-time wave propagations without extreme grid refine- 
ments and small time steps. An integrated approach for the control of numer- 
ical dissipation in high order schemes with incremental studies was initiated 
in [15, 16, 9] and summarized in [17]. Here we further refine the analysis on, 
and improve the understanding of the adaptive numerical dissipation control 
strategy. 

Basically, the development of these schemes focuses on high order nondis- 
sipative schemes and takes advantage of the progress that has been made for 
the last 30 years in numerical methods for conservation laws, such as tech- 
niques for imposing boundary conditions, techniques for stability at shock 
waves, and techniques for stable and accurate long-time integration. We con- 
centrate on high order centered spatial discretizations and a fourth-order 
Runge-Kutta temporal discretizations as the base scheme. Near the bound- 
aries, the base scheme has stable boundary difference operators [12]. To fur- 
ther enhance stability, the split form of the inviscid flux derivatives [5] is 
frequently used for smooth flow problems. To enhance nonlinear stability, 
linear high order numerical dissipations are employed away from disconti- 
nuities, and nonlinear filters are employed after each time step in order to 
suppress spurious oscillations near discontinuities to minimize the smearing 
of turbulent fluctuations. 

Although these schemes are built from many components, each of which 
is well-known, it is not entirely obvious how the different components be best 
connected. For example, the nonlinear filter could instead have been built 
into the spatial discretization, so that it would have been activated at each 
stage in the Runge-Kutta time stepping. We could think of a mechanism that 
activates the split form of the equations only at some parts of the domain. 
Another issue is how to define good sensors for determining in which parts 
of the computational domain a certain feature should be Altered by the ap- 
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propriate numerical dissipation. For the present study we employ a wavelet 
technique introduced in [9] as sensors. Here, the method is briefly described 
with selected numerical experiments. 


2 Base Schemes and Linear Numerical Dissipations 

Consider a conservation law 

n + f{u)x - 0 . 

By a semi-discrete approximation 

^ + Df{uj) = 0, 
at 

where D is a centered finite-difFerence operator, and Uj(t) is the approxima- 
tion of A sixth-order spatial central operator 

Duj = (-Wj-3 + 9uj-2 - 45uj-i -I- 45uj+i - 9uj+2 + WjH-3)/(604\a:) (1) 

is used in our numerical examples. The resulting base scheme with a fourth- 
order Runge-Kutta in time is denoted by CEN66-RK4. Studies also weie 
perform on the fourth and eighth-order centered schemes. The former is less 
accurate and the latter exhibits minor improvement and with a wider stencil 
and more complicated boundary operators over their sixth-order counterpart. 


2.1 Boundary Operators 

On the domain j = 1,2,..., N, the operator (1) traditionally applied at 
the points j = 4,5,...,N - 3 with lower order centered and/or one-sided 
operators for j = 2 , 3 and i = iV - 2. iV - 1. In order to have a stable initial- 
boundary value problem, non-traditional boundary operators that have the 
summation by parts (SBP) property are used. This means that in a weighted 
discrete scalar product, 

{u,v)t,= a,jUiVjAx 

t=l,;=l 

with a positive definite matrix and h = Ax the fixed grid spacing, we 
have the property 

{u,Dv)h = ~{Du,v)h +UNVN - UiVx. (2) 

This implies that whenever stability can be obtained for the PDE by use of 
integration by parts, the same stability estimate can be made for the differ- 
ence approximation by use of (2). Property 2 involved boundary operators 
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for grid points further away from the boundary than the traditional approach 

[ 12 ]. 

In the computations below we use boundary operators having the SBP 
property in a diagonal scalar product, and which are third-order accurate. 
Development of SBP operators was done by Kreiss in the 1970s. Later in [7], 
[12], the theory was revived and extended. 


2.2 Entropy Splitting of the Inviscid Flux Derivatives 

If the above conservation law can be symmetrized by a change of variables, 
u = u(w), then 

u(w)t -f- fw'^x — 0 , 

with = fuUyj symmetric. If we furthermore assume that the flux function 
is homogeneous f(Xw) = A^/(m), with p ^ -1. one can show that Pf = /„m. 
The thermally perfect gas dynamics have a family of variable transformations 
li = u{w,P), defined by u; = which give a homogeneous flux function. 
Here E{u, P) is an entropy function of the conservation law. 

Under the above assumption {p < or /3 > 0 for a perfect gas), the 
flux derivative can be split, resulting in the following form 


Ut + 


/3 + 1 


fx + 


/J + 1 


fwWx = 0 . 


Taking the inner product with w, gives 

-(P+ l)(w,Ut) = P{lU,fx) + (wJwWx) = P{w,fx) + ifwW,Wx)- 
Integration by parts in space and the homogeneity property Pf = fwW, give 
{P T = [m fw'^]a' 


Homogeneity of the change of variables gives 

— {u),u^vj) = (/3 + l)(u(,iu) = —\w fw'^]a- 
Remark: The estimate (3) can be shown to be identical to 

£j'E(u)dx=[Ft 

obtained by integrating the entropy equation, Et Fx — 0, in space. Here F 
is the so called the entropy flux function. 

It follows that wu^w = see [13]. In [5] the entropy splitting is 

extended to conservation laws with nonhomogeneous flux functions. Formulas 
for symmetrizing the compressible Euler equations are given in [3]. Formulas 
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propriate numerical dissipation. For the present study we employ a wavelet 
technique introduced in [9] as sensors. Here, the method is briefly described 
with selected numerical experiments. 


2 Base Schemes and Linear Numerical Dissipations 

Consider a conservation law 


+ f{u)x - 0 . 


By a semi-discrete approximation 

+ £,/(„,) = 0 , 

(it 

where D is a centered finite-difference operator, and Uj{t) is the approxima- 
tion of u{xj,t). A sixth-order spatial central operator 

Duj = (—Uj-3 -t- 9uj-2 — 45uj_i + 45uj+i — 9itj+2 + (1) 

is used in our numerical examples. The resulting base scheme with a fourth- 
order Runge-Kutta in time is denoted by CEN66-RK4. Studies also were 
perform on the fourth and eighth-order centered schemes. The former is less 
accurate and the latter exhibits minor improvement and with a wider stencil 
and more complicated boundary operators over their sixth-order counterpart. 

2.1 Boundary Operators 

On the domain j = the operator (1) traditionally applied at 

the points j = 4,5, ... - 3 with lower order centered and/or one-sided 

operators for j = 2, 3 and j iV - 2, - 1. In order to have a stable initial- 

boundary value problem, non-traditional boundary operators that have the 
summation by parts (SBP) property are used. This means that in a weighted 
discrete scalar product, 

[u^v)^= ^ ^ (Xi jUiVjAx 

t=ij=i 

with <jj j a positive definite matrix and h = Ax the fixed grid spacing, we 
have the property 

(u, Dv)k = ~{Du,v)h + — UiVi . (2) 

This implies that whenever stability can be obtained for the PDE by use of 
integration by parts, the same stability estimate can be made for the differ- 
ence approximation by use of (2). Property 2 involved boundary operators 
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for entropy splitting of a perfect gas case are given in [2], and for a thermally 

perfect gas case are given in [16]. ■ roi i 

The above entropy splitting was used in a numerical method m [2J, where 

the semi-discrete approximation was wiitten as 

•tM = _ - ^Uu,)Dwi. (4) 

Here, P is a difference operator, having the summation by parts property. 
Due to the SBP property, the relation 

— (in, = -wjj } w{wn)'Wn + laf / u; ( w^i (^) 

dt 

in the discrete scalar product follows in the same way aa was shown for 
the PDE. In the special case of periodic boundary conditions, the boundary 
terms disappear, and we conclude that the integrated entropy (w,u^w)^ is 
constant. 

2.3 High Order Linear Numerical Dissipation/Filter 

Numerical dissipation is normally needed in conjunction with entropy split- 
ting. To minimize the lost of accuracy, we employ linear high order numerical 
dissipation obtained by adding an operator of the form 

to the base scheme of order g -f 2 (or as a linear filter after the completion 
of the full step of the time discretization). Here d is a tunable parameter, 
and Ax is the grid spacing. The divided difference operators are defined by 
D+Ui = (uj+i - Uj)/Ax, and = D+Uj-i- The dissipation operator 

has an order of accuracy 2g - 1, and dissipation of order 2g. It is well known 
[6], [8], how to make a boundary modification to obtain an operator which is 
semi-bounded, i,e*, which has the property 

{-iy-\{D+D-Yuj,Uj)h < 0. (6) 

If an extra stabilizing mechanism is desired for difference approximations on 
the split form, it is natural to define the dissipation on the entropy variables. 
The approximation (4) then becomes 

= (L-Df{uj)- -r^fw(uj)Dwj+d{-l)‘‘~'^Ax^‘’''^{D+D-)Hvj. 

dt /3+1 '' '' i9 + l 

If the SBP property and the semi-boundedness for the dissipation both hold 
in the same scalar product, we take the scalar product with Wj, and obtain 
as previously 
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^{w,Uw'^)h. = — ~ M\h- (^) 

at 

The dissipation thus helps to control the increase of the entropy. An example 
where this applies is for periodic problems. For general boundary conditions, 
the dissipation has to be modified somewhat in order to obtain (8), because 
the scalar product for the summation by parts property is not the same as 
the standard scalar product used for the dissipation estimate. For the sixth- 
order base scheme, eighth-order linear dissipation {ci = 4) is used and is 
denoted by CEN66-D8 without entropy splitting, and CEN66-ENT-D8 with 
the splitting. 


Isentropic Vortex Evolution 

{Honzontally Convecting Vortex, vortex strength ^=5} 
Freestream : 

(uoo.Wo.) = (1,0); p»- ftm, = I 

1C .' Pcflurtalloni are added lo the IrMiUtam (not in enlrap)^) 

A , 

(Su,Sv] = _ ,.)) 

ST = , ( . 7~ 

Sjir* 

Computational Domain t Grid Size : 

0<a<lfllt-5<v<5 
80 X 79 Unitonn grid 
Periodic BC in z & y 

Fig. 1. Vortex convection problem description. 
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3 Long-Time Integration of Smooth Flows 

Extensive numerical experiments on a 2-D compressible inviscid vortex con- 
vection with or without entropy splitting and/or nonlinear filter were con- 
ducted in [16, 9]. They showed that entropy splitting alone and/or nonlinear 
filter helps stabilize the scheme for a much longer convection time with al- 
most perfect vortex preservation. Here we study the same problem for the 
behavior of the high order linear numerical dissipation and compare the so- 
lution with other well known methods. This is a pure convection problem 
and the exact solution should retain its shape. The problem is described in 
Fig. 1. The boundary conditions are periodic, and the solution is smooth, so 
that difficulties coming from boundary conditions and from non-smoothness 
of the solution are avoided. The computational domain is 10 x 10 with a 
uniform 80 x 80 grid. The time scale is such that one period corresponds to 
10 time units. Density contours without entropy splitting and without any 
added dissipation (CEN66-RK4), are shown in Fig. 2. The solution after 5 
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Fig. 2. Density contours after 0, 5, and 5.6 periods using CEN66-RK4. 


periods has developed small, hardly visible, oscillations, and the computation 
breaks down before reaching 6 periods. The oscillations are due to nonlinear 
instabilities. Density contours using the entropy split form (4) (CEN44-ENT- 



Fig. 3. Density contours using CEN66-ENT-RK4 with entropy split parameter, 
/?=!. 

RK4) is shown in Fig. 3 with j3 = 1 (equal weight on the conservative and 
the non-conservative part). For this problem, P = 1 appears to be the best 
choice. Almost perfect vortex preservation is nearly 10 times longer than the 
unsplit case. However, the solution eventually breaks down (after 68 peri- 
ods). For the entropy split method, it follows from (5) that the integrated 
entropy should be constant. Figure 4 shows the entropy integral as function 
of time for the computation in Fig. 3. The integrated entropy does decrease. 
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and the estimate holds, but is too weak. However, we have found that there 
are other values of /3 where the computed entropy integral increases. 0 does 
not appear in the semi-discrete estimate, but is apparently important for the 
time-discrete computations. It can be seen that entropy splitting retain the 


Norm vs time 



Fig. 4. Entropy integral vs. time of CEN66-ENT-RK4 with d — ???, 0 = 


vortex shape from 5 periods to approximately 40 periods. However, in many 
applications, for example computation of rotorcraft flows, one needs to follow 
vortices for several hundred periods. We next investigate whether the use of 
high order linear numerical dissipation can be used to further increase the 
number of periods. Figure 5 shows the density contours of the sixth-order 



Fig. 5. Density contours by CEN66-D8-RK4 after 30, 150 and 200 periods. 


base scheme with the eighth-order numerical dissipation added (CEN66-D8- 
RK4). The solution is very accurate up to approximately 150 periods, but 
the accuracy degenerates at later times. Although the accuracy becomes very 
poor, the solution did not break down. The computation could be run up to 
300 periods, the maximum time we used. The coefficient d in front of the 
dissipation operator influences the behavior of the solution. The errors of 
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the solution as function of time for different values of the dissipation param- 
eter d are shown in Fig. 6. The error norm is computed as the sum of the 
error norms of the four components. From Fig. 6, we conclude that there is an 
initial phase where the lower the dissipation, the better the accuracy, and the 
error increases in a “regular” way. At later times the larger dissipations have 
two points where the error increases quickly. For example, for d = 8 x 10 
the first error increase comes at 70 periods and the second at 190 periods. 
We also conclude that for long-time integration, using the smallest possible 
dissipation is not optimal. After 120 periods the dissipation 1 x 10 ^ gives 
a more accurate solution than the dissipation 5 x 10 The smaller dissipa- 
tion is insufficient to suppress the nonlinear instabilities. After 200 periods, 
all dissipations end up in a state where all accuracy is lost. The norm of 
the difference between the solution and the constant state with no vortex is 
approximately 0.1. Results using both entropy splitting and the linear nu- 
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Fig. 6. Error vs. time for different dissipation coefficient d for CEN66-D8-RK4. 


merical dissipation CEN66-ENT-D8 exhibit very minor improvement and no 
improvement at larger period of the vortex convection over the CEN66-D8 
case. A possible explanation is that the very delicate cancellation properties, 
which keep the entropy split method stable, are destroyed by the numerical 
dissipation as time progresses. A fifth order WENO scheme (WEN05-RK4) 
[4], the so called DRP scheme [14], and a second-order accurate MUSCL 
TVD scheme (MUSCL-RK2) were computed for comparison. The WEN05 
scheme is a general purpose scheme suited for flows with shock waves with 
no parameters to tune. The DRP scheme is designed for aeroacoustics. It is 
fourth-order accurate in space and second-order in time. Nonlinear instabili- 
ties are suppressed by a low order numerical dissipation term with a tunable 
parameter. We have tried to choose this dissipation strength as small as pos- 
sible while keeping the scheme stable. Density contours are shown in Fig. 
7 after 150 periods, and in Fig. 8 after 200 periods using the same contour 
levels. The MUSCL-RK2, and the DRP scheme are extremely diffusive. Both 
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(a)IC 




(b) MUSCL 
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(c) DRP 





(dpEN66-D8 



(e) WEN05 


Fig. 7. Comparison of different methods: Density contours at 150 Periods 



(a) IC 



(b) MUSCL 



(c) DRP 




(e) WEN05 


Fig. 8. Comparison of different methods: Density contours at 200 Periods 


the WEN05 and the CEN66-D8 give better results. The error as function of 
time for the different methods is shown in Fig. 9. The error of MUSCL is not 
shown, instead the dashed curve shows the error of the entropy split sixth- 
order scheme, with /3 = 1, and no numerical dissipation (CEN66-ENT). After 



Periods 


Fig. 9. Error vs. time for different methods. 
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200 periods, the result from the WEN05 scheme looks more pleasing to the 
eye, but as shown in Fig, 9, its error norm is larger than that of the centered 
scheme. The centered scheme loses accuracy due to dispersive errors, whde 
for the other schemes, diffusive errors dominate. The performance of nonlin- 
ear filter scheme (including the use of the dissipative portion of WEN05 as a 
nonlinear filter) computations (same base scheme) are slightly more diffusive 
than CEN66-D8, and is reported in [9, 16, ?]. Studies on the blending of the 
linear dissipations with the nonlinear filters are reported in [17]. 


4 Nonlinear Filter/Numerical Dissipation for 
Discontinuities 


The linear numerical dissipation used in the previous section is not strong 
enough to suppress the spurious oscillations due to the presence of shock and 
shear waves. A filter approach using nonlinear dissipation as a post processing 
filter, which is applied after each time step, is employed [15]. The filter consists 
of a conservative step 

where u* is the solution computed by the high order base scheme at time level 
n-(- 1, and A is the ratio of time step over Ax. For the numerical example the 
filter flux h*is obtained from a second-order TVD scheme by subtracting the 
centered difference part and inserting a sensor, i.e., 

^j + l/2 - ~ (/j + 1 + 

The flux numerical flux of a second-order accurate TVD scheme. 

This TVD^-flux typically contains flux limiters and Riemann solvers. The 
sensor s,+i /2 is made such that the dissipation is switched on only near 
discontinuities. One example is to compute the local wavelet coefficients of 
the flow field by letting s,+i /2 be close to unity where the wavelet coefficients 
show that the solution has poor regularity, and taking Sj+ 1/2 to be zero 
otherwise. See [9] for details. Computations using the dissipative portion 
of the WEN05 schemes or seventh-order high-resolution schemes were also 
considered. See [18] for an example using the WEN05 nonlinear dissipation. 


5 2-D Shock/Shear /Boundary-Layer Interaction 

Extensive grid refinement for shock-turbulence interactions, a complex multi- 
scale viscous combustion flows and a shock/shear/boundary-layer interaction 
are reported [17, 10, 11, 1]. Here we illustrate selected results for the complex 
shock/shear/boundary-layer interaction. An ideal gas is at rest in a 2-D box 
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0 < x^y < 1. A planar shock of Mach 2.37 located at x — 1/2 at time zero 
separateslwo different states of the ga5. The 2-D compressible Navier-Stokes 
equation with no slip BCs at the adiabatic walls is used. The solution will 
develop complex 2-D shock/shear/boundary-layer interactions, which depend 
on the Reynolds number. The dimensionless initial states given in [1] are 

Pi = 120, pz, = 120/7, Pfi = 1-2, Pfi = 1 - 2 / 7 i 

where pi,pL are the density and pressure respectively to the left of x = 1/2, 
and PR, PR are the same quantities to the right of x = 1/2. All velocities are 
equal to zero, 7 = 1.4, the Prandtl number is 0.73, and the Reynolds number 
is 1000. All the computations stop at the dimensionless time 1 after the shock 
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(a) 1000 X 500 grid 
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(c) 3000 X 1500 grid 


(d) 4000 X 2000 grid 


Fig. 10. Grid refinement study: Density contours of MUSCL RK2. 


wave has reflected from the right wall and has almost reached the middle of 
the domain, x = 1/2. Due to the rapidly developing symmetry flow structure, 
a uniform Cartesian grid with the lower half of the domain is used. A grid 
convergence study of the sixth-order method with wavelet as the sensor to 
filter the solution using the nonlinear numerical dissipation (WAV66-RK4) is 
shown in Fig. 11. The same computation using the MUSCL- RK2 is shown in 
Fig. 10. Except for the largest vortical structure, the 1000 x 500 grid solution 
of WAV66-RK4 captures the overall flow structure of a 4000 x 2000 grid. 
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whereas the MUSCL scheme is still not convergent ever at a finer grid of 
4000 X 2000. Since the filter is only applied once after each time step, the 
computational cost of WAV66-RK4 is almost the same as the computational 
cost of MUSCL-RK2 with the added advantage that the physical viscosities 
were taken into consideration by WAV66-RK4. 


(c) 3000 X 1500 grid (d) 4000 x 2000 grid 




WAvae 1000x500 



(a) 1000 X 500 grid 


WAV66 2000x1000 



(b) 2000 X 1000 grid 


Fig. 11. Grid refinement study: Density contours of WAV66-RK4. 
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